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In cells, most of cargos are transported by motor proteins along microtubule. Biophysically, 
unidirectional motion of large number of motor proteins along a single track can be described by 
totally asymmetric simple exclusion process (TASEP). From which many meaningful properties, 
such as the appearance of domain wall (defined as the borderline of high density and low density of 
motor protein along motion track) and boundary layers, can be obtained. However, it is biologically 
obvious that a single track may be occupied by different motor species. So previous studies based 
on TASEP of one particle species are not reasonable enough to find more detailed properties of 
the motion of motors along a single track. To address this problem, TASEP with two particle 
species is discussed in this study. Theoretical methods to get densities of each particle species are 
provided. Using these methods, phase transition related properties of particle densities are obtained. 

Our analysis show that domain wall and boundary layer of single species densities always appear 
simultaneously with those of the total particle density. The height of domain wall of total particle 
density is equal to the summation of those of single species. Phase diagrams for typical model 
parameters are also presented. The methods presented in this study can be generalized to analyze 
TASEP with more particle species. 
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Many driven diffusion systems have been developed 
to model intracellular motility [1, 2]. Among which, 
one-dimensional totally asymmetric simple exclusion pro¬ 
cess (TASEP) is usually employed to describe the uni¬ 
directional motion of large number of motor proteins 
along microtubule, the transcription process of RNA 
polymerases along DNA, and the translation process of 
ribosomes along messenger RNA transcript [3-5]. In 
TASEP, motion track of particles is simplified to be a 
one-dimensional lattice with length TV + 1, particles en¬ 
ter the track at first/initiation site 0 provided site 0 is not 
occupied, and leave track from last/termination site TV. 
If lattice site i + 1 is not occupied, the particle at lattice 
site i will hop forward to site i +1 with given rate. Gener¬ 
ally, particles may also detach into environment and end 
their unidirectional motion from bulk lattice sites, and 
new particles may attach to any one of unoccupied bulk 
sites. In last decades, TASEP has been studied exten¬ 
sively, especially for phase transition related properties, 
i.e., the appearance of domain wall (DW) and bound¬ 
ary layer (BL), which are usually driven by boundary 
conditions [6-12], However, in most of previous stud¬ 
ies, the self-propelled particles are usually assumed to be 
from the same species. But in cells, one protofilament 
of microtubule may be occupied by motor proteins from 
different species [13]. So, to know more detailed proper¬ 
ties about the motion of motor proteins in real cells, the 
usual TASEP should be generalized to include particles 
from different species [14-17]. 
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The simplest generalization is two-species ASEP, 
which has been discussed in recent studies [18-21]. How¬ 
ever, in almost all previous studies, no particle attach¬ 
ment/detachment is allowed to/from bulk sites of the 
track, i.e. the site T for 1 < i < TV — 1. They usu¬ 
ally assumed that one species enters track at site 0 and 
leaves from site TV, while the other species enters at site TV 
and leaves from site 0. Meanwhile, it is also assumed that 
two-species pair (Pi)j(P 2 )i+i can change to (P 2 )i(Pi)i+i, 
i.e. forward hopping of one species is not blocked by the 
other species. Here (Pk)i means there is a particle Pk 
at site i. For convenience, the two particle species are 
denoted by Pi and P 2 respectively. 

The TASEP discussed in this study also includes two 
particle species. But both of them enter track at ini¬ 
tiation site 0 and leave from termination site TV, i.e. 
they not only travel along the same track, but also 
move to the same direction. Previous studies about 
one species TASEP have shown that nontrivial attach¬ 
ment/detachment is one of the key driven factors to 
the appearance of DW in particle density along track 
[4, 7, 11, 22-24]. So, in this study, both of the two species 
are allowed to attach to (and detach from) bulk sites of 
the track. The same as in one species cases, we also called 
this process “TASEP-LK” process [7]. 

Let Hi and rrii be occupation numbers of species Pi and 
P 2 at site i, respectively. Specifically n, = 1 means site 
i is occupied by a particle Pi, while rii = 0 means site i 
is not occupied by particle Pi. Because of the hard-core 
exclusion, n, + rrii = 0 or 1. For 1 < * < TV — 1, the 
time evolution of rii and rrii are governed by following 
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equations 

drii/dt = n*-i(l - n» — m») - nj(l - n»+i - m»+i) 

+£*> i , a (1 - ni — mi) - uJi tD ni, ( 1 ) 

drrii/dt = qmi-ffl - ni — mi) — gm»(l — n i+ i — m i+ i) 

+w 2 ,a( 1 - «» - r/ij) - U2 }D mi. (2) 

Where, for convenience, forward hopping rate of species 
Pi is normalized to be unit, and forward hopping rate of 
species P 2 is assumed to be q < 1. The attachment rate 
of species P k to any unoccupied bulk site i is denoted by 
uik } Ai and the detachment rate of species Pk from bulk 
sites is denoted by u) k ,D- At initiation site i = 0 (left 
boundary), 

dno/dt = ai(l — no — mo) — no(l — ni — mi), (3) 
dmo/dt = a 2 (l — no — mo) — qmo(l — n\ — mi). (4) 

Where ak is the entry rate of particle Pk from environ¬ 
ment. While at termination site i = N (right boundary), 

dn^/dt = njv—i(l — tin — mjv) — PiUn, (5) 
dm^/dt = qrriN-i(l — tin — mw) — fcmN, (6) 

with [dk the leaving rate of particle Pk to environment. 

In cells, hopping rate of particles is usually deter¬ 
mined by their biochemical properties. For example, for¬ 
ward hopping of motor proteins, such as conventional 
kinesin, is mechanochemically coupled with ATP hydrol¬ 
ysis. Each forward mechanical step is tightly coupled 
with one ATP hydrolysis. Therefore, hopping rate is 
determined by the rate of ATP hydrolysis. Meanwhile, 
leaving rate of particles from motion track is also deter¬ 
mined by their biochemical properties, or even the rate 
of ATP hydrolysis. Therefore it is biophysically reason¬ 
able to assume that the ratio of leaving rate of the two 
species /3i : /3 2 is equal to ratio of their forward hop¬ 
ping rate 1 : q. Meanwhile, experiments found that, if 
the two heads of motor protein kinesin are both in ADP 
binding state, then it will soon detach from microtubule. 
This means that detachment rate of motor proteins are 
also determined by their biochemical properties and the 
rate of ATP hydrolysis. So, for convenience of theoretical 
analysis, this study also assumes that u>\ . u '■ w 2 ,d = 1 : q- 
However, except their biochemical properties, entry rate 
ak and attachment rate uJk,A of species Pk are also influ¬ 
enced by environmental conditions, especially their con¬ 
centrations. So, the corresponding rate ratios may be 
different from the ratio 1 : q of forward hopping rate. 

Defining p\{i) = ( ni) and p 2 (i) = {mi). By mean-field 
approximation, for large track length N limit, equations 
for steady state values of densities p\ and p 2 can be ob¬ 
tained from Eqs. (1,2), 

d x Jk = g 1-fc f4,A(l - Pi ~ Pi) ~ ttp k , k = 1,2. (7) 

Where J k = p fe (l-pi-p 2 ), ^ k,A = Nu k ,A, ^ = Nui, d = 
A r tU 2 ,D/q ='■ N/uj, and 0 < x < 1. Note, this study 
assumes that (nimi+i) = p\{i)p 2 {i + 1), and in MFA 


analysis the track length is normalized to be 1. From Eq. 
(7), one can show that total particle density p = pi + P 2 
satisfies 

d x J = KQ{ 1 - p) - ftp, (8) 

where K = K\ + AT 2 = [u\ z a + ^ 2 ,A/q]/^, and J = p{ 1 — 
p). From Eqs. (3-6), one can show that at boundaries 
x = 0,1, particle densities pi,p 2 and p satisfy pi(0) = 
aq,/O 2 ( 0 ) = a 2 /q, and p(l) = pi(l) + p 2 { 1) = 1-/3 with 
/3 = /3i = P 2 /q- 

Eq. (8) implies that the governing equation for total 
density p is the same as the one in one-species “TASEP- 
LK” process [26]. But with effective detachment rate 
fljD = fl, effective attachment rate £Ia = KQ, 1 initiation 
(entry) rate a = a± + a 2 /q, and termination (leaving) 
rate /3 = /3i = /3 2 /q (see Fig.l). So total density p can be 
obtained by the same method as in one-species “TASEP- 
LK” process. But the main difficulty for two-species cases 
is how to get single species densities pi and p 2 ■ Actually, 
no reasonable boundary conditions at x = 1, i.e., values 
of p\{l) and p 2 (l), can be derived from Eqs. (5,6). Mean¬ 
while, properties of densities p\ and p 2 are different from 
the ones in one-species cases. For example, in one-species 
“TASEP-LK” process, particle densities before and after 
DW location x w satisfies p{x~) + p{x+) = 1. This is 
because that across location x w , current J = p{ 1 — p) is 
conserved. But for two-species process, particle densities 
Pi, p 2 do not satisfy this relation. Instead, conservation 
of current J k = p*, (1 — p) gives that across DW location, 
density ratio p\/p 2 is not changed (see Fig.l). The plots 
in Fig.l imply that DWs of density p,pi,p 2 appear at 
the same location. Further numerical calculations show 
that their BLs also appear simultaneously, see Fig. S4 in 

[25]. 



FIG. 1: (Color online) Mean density profiles obtained from 
Eqs. (l-6)(black dotted line), and their mean-field approxi¬ 
mations obtained from Eqs. (7,8) (cyan solid lines for p, red 
dashdot lines for pi, and blue dashed line for p 2 ). Parameter 
values used in calculations are (a) Di,a = 0.01, IH.a = 0.14, 
an = 0.2, a 2 = 0.09, and (b) Di,a = 0.1, 0,2,a = 0.05, 
ai = 0.1, o 2 = 0.18. Other parameter values are N = 10 3 , 
q = 0.9, ni,D = O 2 ,o/q = 0.1, /3i = /3 2 /q = 0.3. Total par¬ 
ticle density p in both (a) and (b) is the same as the one in 
usual “TASEP-LK” process with a = /3 = 0.3, and Od = 0.1, 
Oa = 0.15. Locations of domain wall for total density p and 
single species densities p\,p 2 are the same. Across domain 
wall location, current J = p(l —p) and current J k = pfc(l — p) 
for k — 1,2 are all conserved. 
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Along motion track, density currents J, J\ , J 2 al¬ 
ways change continuously, even if corresponding den¬ 
sities p, pi,p 2 are discontinuous. At any location x 
of the track, p(x~)[ 1 — p{x~)\ = J{x~) = J(x + ) = 
p{x + )[l - p{x + )] 1 and p k {x~)[l - p(x~)\ = J k {x~) = 
J k {x + ) = p k (x + )[ 1 — p{x + )]. Through simple anal¬ 
ysis, we obtained p k {x~)/p{x~) = p k (x + )/p(x + ) or 
Pi{x~)/p2{x~) = pi{x + )/p2(x + ). This relation also 
holds at DW location x w . 

As we have mentioned before, total density p can 
be obtained from Eq. (8) with boundary conditions 
p(0) = pi(0) + p 2 (0) = ai + a 2 /q and p( 1) = 1-/3. 
But, if DW exists between boundaries x = 0 and x = 1, 
then without boundary value p k { 1), density p k cannot be 
directly determined by Eq. (7). Actually, with bound¬ 
ary value p k ( 0), only the value of density p k before DW 
location x w can be directly obtained by Eq. (7). One of 
the main aims of this study to find methods to get single 
species density p k along the whole track from Eq. (7), 
but with only boundary values at x = 0. 

Let A := p(x+) — p( x~), i.e. the DW height 
of total density p. Then A k := p k {x+) — 
Pk(x~) = [Pk{x+)/p{x+)}p(x+)-{p k (x~)/p{x~)}p(x~) = 
[Pk{x+) / p(x+)\A = [p k (x-)/p(x~)]A. So, A = Ai+A 2 , 
and A > 0 iff A k > 0. Therefore, DWs of total density 
p and single species densities p\ and P2 always appear at 
the same location, which is consistent with the finding in 
numerical calculations (see Fig. 1). 

The value p k (x+) can be obtained as follows, p k {x+) = 
Pk(x~) + At = Pk(Xw)+[pk(x-)/p(x-)][p(x+)-p(x-)] = 
P( x t,)Pk{x w )/p(x w ). Using p k {x+) as left boundary con¬ 
dition, the single species density p k after DW location x w 
can be obtained from Eq. (7). Note, the above method 
to get single species density p k is also applicable for the 
cases where there exists BL in density p. BLs of single 
species densities pi, P2 and total density p are also appear 
simultaneously. 

If K = 1, Eq. (8) reduces to 

(d x p-Q)(2p-l) = 0, (9) 

which yields two solutions. The constant solution p = 
1/2 coincides with the density pi = K/(K + 1) given by 
Langmuir kinetics (LK), and corresponds to the maximal 
current (MC). The other solution which matches left or 
right boundary condition is p a (x) = fix + a or pp (x) = 
flx+l — ft — fl. Using these three solution candidates and 
based on the continuity of currents J, Ji, J 2 , total density 
p and single species densities p \, P2 can be obtained [25] . 
Typical examples of p, pi, p 2 for K = 1 are plotted in Fig. 
2, including the cases of low density (LD) phase {p < 
1/2), high density (HD) phase (p > 1/2), the maximal 
current (MC) phase (p = 1/2), as well as left and right 
BLs. Where thick dotted lines are obtained by numerical 
iterations of Eqs. (1-6), and others are obtained from 
Eqs. (7,8) using the method presented here [25]. Fig. 2 
shows that for large initiation rate a = ot\ + a 2 > 0.5, 
left BL exists in both total density p and single species 
densities pi and P2 ■ Meanwhile, right BL may appear 


for either large or small values of termination rate /3. 
Again, DW of single species density p k occurs at the same 
location as that of total density p. Which satisfies A k = 
[p k (x-)/p(x-)\A and A = Ai + A 2 . 



X 

FIG. 2: (Color online) Examples of density profiles for K = 1. 
The same line types as in Fig. 1 are used, (a) LD-BL phase 
with a < 0.5 and /3 < 0.5; (b) BL-HD phase with a > 0.5 and 
/? < 0.5; (c) LD-DM-HD phase with a < 0.5 and /3 < 0.5; 
(d) LD-MC-BL phase with a < 0.5 and /3 > 0.5; (e) BL- 
MC-HD phase with a > 0.5 and /3 < 0.5; (f) BL-MC-BL 
phase with a > 0.5 and /3 > 0.5. Where a = on + 02 /( 7 , 
and P = Pi = p 2 /q- Unless explicitly presented, in all the 
following calculations q = 0.9 is used. 

For special cases K = 1, phases of single species den¬ 
sities pi and p 2 may not be the same as that of total 
density p. If p is in LD phase, then both pi and p 2 will 
be in LD phase, see Fig. 2(a). However, when p is in HD 
phase, one of the single species densities, or even both of 
them, may still remain in LD phase, see Fig. 2(b,e). For 
special cases K = 1, if initiation rate a and termination 
rate /3 satisfy [25]. 

a > 1/2 —fl, /3 > 1/2 —D, and a+/3 > 1-D, (10) 

then density p may be in maximal current phase (MC) 
near right boundary x = 1, in which p = 1/2 and cur¬ 
rent J = J\ + J 2 = 1/4. Depends on rates a k and fl k ,A , 
single species density p k may not be constant, and they 
may increase or decrease along motion track. This is dif¬ 
ferent from the cases in which density p is in LD or HD 
phase. For those cases both pi and p 2 increase with x , 
and with slops and fl2,A/q, respectively, see Eqs. 
(S7,S7) in [25]. It can be shown that when p = 1/2, 
slopes of pi and p 2 have same absolute value, but op¬ 
posite signs. If Di^o^ — fl 2 ,AOti > 0 then pi has pos¬ 
itive slope, i.e. increases along the motion track, while 
p 2 has negative slope. Both pi,p 2 are also constants iff 
fli,A &2 ~ fl 2 ,AOti = 0, see Fig. S2 in [25]. 

Based on current continuity, and using expressions of 
p a and pp, which satisfy Eq. (7) with boundary condi¬ 
tions Pa (0) = a and pp(l) = /3 respectively, we found 
that for AT = 1, DW appears iff \a — /3| < fl and 
a+P+fl < 1 [25], The DW lies at x w = (D+/3-a)/(2D) 
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with height A = 1 — a — ft — fl. The DW height of density 
Pk can then be obtained by A^ = [pk{x~) / p(x~)\ A. see 
Eqs. (S4,S5,S7,S8) in [25], 

Using similar method as in [26], typical examples of 
phase diagram of density p, for special cases K = 1, are 
given in Fig. 3. Generally, density p may be in any one of 
the seven possible phases, which include all combinations 
of LD, MC, and HD. The lines in Figs. 3(b,c,e,f) are 
obtained by 0:2 = q(a — a 1 ), since density p depends only 
on the summation a = a\ + a^/q- 



FIG. 3: (Color online) Phase diagrams of total density p = 
pi + pi obtained from stationary mean-field Eq.(8) with K = 
1, Q = 0.3 in (a,b,c), and 12 = 0.5 in (d,e,f). (b,c) and 
(e,f) are phase diagrams in ( 01 , 02 ) plane, which correspond 
to ft = 0.4 and ft = 0.6, respectively (see the horizontal dotted 
and dashed lines in (a,d)). 

As mentioned in [7], due to particle-hole symmetry we 
only need to discuss the “TASEP-LK” process for AT > 1. 
For K > 1, the idea used in the special cases K = 1 can 
also be employed to get single special density pk [25]. 
Roughly speaking, total density p, as well as its loca¬ 
tions of DW and BL, can be obtained from Eq. (8) 
with boundary conditions p(0) = a = a± + ai/q and 
p(l) = 1 — ft = \ — ft\ = 1 — /? 2 /9- Then single species 
density pk in interval [0,Xu,) can be obtain from Eq. (7) 
with left boundary condition aq or 0 : 2 /q- Finally, den¬ 
sity pk in interval ( x w , 1] can be obtained from Eq. (7) 
with left boundary condition pfc(xj), which is given by 
Pk{Xy,) = Pk{x~) + Afc. The main difference from that in 
special cases K = 1 is that Lambert function [27] should 
be employed to help to get densities p and pk [25]. 

For K > 1, properties of total density p are similar 
as the ones of the one-species “TASEP-LK” process [7, 
11 ], i.e. there may exist left or right boundary layer, 
domain wall, or “Meissner” (M) phase, see Fig. 4. Here 
“Meissner” phase means that density p satisfies 1/2 < 
p < pi = K/(K +1), and is independent of initiation 
rate a and termination rate ft, see Fig. 4(f). Properties 
of single species density pk may be different from p. For 
example, if p lies in [0,1/2] or [pi, 1], it will increase along 
the track. Otherwise, p decreases along the track. This is 
because that dp = (AT+l)D(p— pi)/(2p— 1), see Eq. (S20) 


in [25]. Therefore, if there exists DW in (0,1) and the 
termination rate 1 — pi < ft < 1, then after DW location 
x w , density p decreases monotonically. Otherwise, if 0 < 
ft < 1 — pi, p will increase after x w . However, the results 
in Fig. 4(b,c) show that, after DW location x w , the 
monotonicity of pk may be different from p. 
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FIG. 4: (Color online) Typical examples of density p for 12 = 
0.1 and K = 1.5, i.e. pi = K/(K +1) = 0.6. (a) LD-DW-HD 
phase, with a < 0.5 and ft > 0.5; (b) LD-DW-HD phase, with 
a < 0.5 and 1 — pi < (3 < 0.5; (c) LD-DW-HD phase, with 
a < 0.5 and 0 < ft < 1 — pr, (d) HD phase, with a < 0.5 and 
ft < 1 — pi; (e) LD phase, with a < 0.5 and 1 — pi < ft < 0.5; 
(f) “Meissner” phase, with a > 0.5 and ft > 0.5. Line types 
are the same as in Fig. 1. 

Previous analysis about special cases K = 1 has shown 
that, if total density p is in MC phase, i.e. p = 1/2, then 
single species density pk may not be constant. Similar 
results hold for the general K > 1 cases. If /3 = 1 — p*, 
then p = pi = K/(K+ 1) is constant near right boundary 
x = 1. The plots in Fig. S5 show that the corresponding 
single species densities pk may not be constant. Theo¬ 
retical analysis gives that the monotonicity of pk , when 
p = pi is constant, is also determined by the sign of 
^■i.ao -2 ~ ^ 2 , 7 !op. With positive values of it, pi increases, 
while P 2 decreases, along motion track. Both densities pi 
and P 2 will be constant iff — ^ 2 ,ao.i = 0, see [25] 

for detailed analysis. 

Examples of phase diagram of total density p in (a, /?) 
plane, for general cases K > 1, are plotted in Fig. 
5(a,b,c). Similar as in the one-species “TASEP-LK”, BL 
may appear at one or both of the two boundaries. DW 
may appear in interval (0,1), and density p may be in LD 
phase (p < 1/2) or HD phase (p > 1/2). To show more 
details about the “TASEP-LK” process, in Fig. 5 the 
HD phase is divided into two different cases, HDi phase 
(1/2 < p < pi = K/{k + 1)) and HD 2 phase (p > pi). 
From the phase diagram in (a, ft) plane, phase diagrams 
in any planes of parameter pair (tri, 172 ) can be easily ob¬ 
tained. Where <Jk = «i, « 2 , fti, /? 2 - Examples of phase 
diagram in (au,^) plane are plotted in Fig. 5(d,e,f), 
which are corresponding to the dotted horizontal lines in 
Fig. 5(a,b,c) respectively. The above discussion about 
the relationship between total density p and single species 
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FIG. 5: (Color online) Phase diagrams of total density p ob¬ 
tained by stationary mean-field Eq.(8) for K > 1. Parameter 
values used in calculations are fl = 0.1, K = 1.5 for (a,d); 
fl = 0.15, K = 1.5 for (b,e); and Q. = 0.1, K = 4 for (c,f). 
(d,e,f) are phase diagrams in (ai, « 2 )-plane with termina¬ 
tion rate f3 = 0.45,0.7,0.3, i.e., corresponding to horizontal 
dotted lines in (a,b,c), respectively. The meanings of Roman 
numerals in each connected area of figures are as follows. (1) 
I:=LD-BL+, which means that density p is in LD phase, and 
there exists right boundary layer (BL) in which p increases 
sharply, i.e. p(l) > p(l — e) for small real number e. (2) 
II:=LD-BL“, i.e., p is in LD phase, and it decreases sharply 
at right boundary. (3) III:=LD-DW-HDi, near left boundary 
density p is in low density phase, while near right boundary p 
is in high density phase, and there exists DW between these 
two phases. Here “HDi” means p > pi, i.e., /3 < 1 — pi. 
(4) IV:=LD-DW-HD 2 , where HD 2 means that density p is 
between 0.5 and pi, i.e., I — pi < /3 < 0.5. (5) V:=LD- 

DW-M r , where M r means that density p is between 0.5 and 
pi but its value is independent of right boundary condition, 
therefore right BL appears. M r can be regarded as the right 
half part of “Meissner” phase. For this phase, /3 > 0.5. (6) 
VLssBL^-HDi, i.e., p > pi and p increases in left BL. (7) 
VII:=BL) f -HD 2 . (8) VIII:=M + , which means that p is in 
“Meissner” phase and it increases in left BL. (9) IX:=BL ( “- 
HDi. (10) X:=BL”-HD 2 . (11) XI:=BL,--M. VIII and XI 
are two species cases of “Meissner” phase. For examples of 
density p which is in one of the above eleven phases, see Fig. 
S4 in [25]. 


density pk implies that phase diagrams of density pk are 
the same as those of the total density p. 

In summary, “TASEP-LK” process with two particle 
species is discussed in this study. Different from previ¬ 
ous studies about two-species TASEP, particle attach¬ 
ment/detachment to/from bulk sites of motion track is 
allowed, both of the two particle species enter into the 
track from the same boundary, and move unidirection- 
ally to the same direction. The two particle species do 
not change to each other, and do not change their lo¬ 
cations even if they are adjacent to each other [18-21]. 
This study found that, domain walls and boundary lay¬ 
ers of total particle density and single species densities 
always appear simultaneously. The height of domain wall 
of total density is equal to the summation of those of the 
two single species. Based on these properties, theoretical 
methods to obtain steady state densities of the two par¬ 
ticle species are presented. Our results show that prop¬ 
erties of single species densities may be different from 
that of the total particle density. In this study, phase 
diagrams of particle density in typical parameter planes 
are also presented. The methods presented in this study 
are also available to the analysis of “TASEP-LK” process 
including more than two particle species. The results of 
this study are helpful to further understandings of the 
biophysical process of cargo transportation in living cells, 
where one single protofilament of microtubule is actually 
occupied by various kinds of cargos and motor proteins 
[13, 28], 
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